'''
File: calculateGeneMetrics.py
Author: Adam Hockenberry
Description: I am accumulating different ways to quantify the codon usage for a single gene here.
The functions towards the end are predominately ways of testing that everything is working but 
the functions towards the beginning all accept a basic codon list and can output a single metric.
All methods assume one standard codon usage table and the codon list should not contain start or stop codons.
'''

#Standard path imports
from optparse import OptionParser
from math import log, exp
from scipy import stats
import numpy as np

import sequenceAnalysisFunctions as saf
import parsingGenbankFiles as pgb
from matplotlib import pyplot as plt
from Bio import SeqIO

#Non-standard imports

import codon_entropy
import geneMetrics
import calculate_CAI
import random
import compileExpressionData
import os
import glob
from scipy.spatial import distance
import operator
import cod_trna_aa_entropy as cod_trna_aa_entropy

import beeswarm
from mpltools import style
from mpltools import layout

style.use('ggplot')
#Global directories and variables

def growth_metric_simple():
    genomeDicty = {}
    f = open('../Data/growthRatesConvertedEdited.csv').readlines()
    for x in f:
        listy = x.split('\t')
        print listy
        listy[3] = listy[3].strip('\n')
        if len(listy[3]) > 5 and len(listy[2]) > 0:
            genomeDicty[listy[3].strip('\'').strip('"')] = float(listy[2])
    metricWright = []
    metricNovembre = []
    metricInfoGC = []
    metricTRNA = []
    metricCodGivenTRNA = []
    genTimes = []
    for x in genomeDicty.keys():
        if x.find('Mycoplasma') == -1:
            print '########################################################################' 
            print x
            print '########################################################################' 
            genes = []
            lociNames = []
            for genbankFile in glob.glob(x + '*.gbk'):
                genesSingle, lociNamesSingle = pgb.get_coding_sequences(genbankFile, returnLociNames=True)
                if len(genesSingle) > len(genes):
                    genes = genesSingle
                    lociNames = lociNamesSingle
            genesDict = {}
            for y in range(len(lociNames)):
                genesDict[lociNames[y]] = genes[y]
            print "I am using", len(genesDict.keys()), "genes for this organism"
            codonList = [saf.split_string_into_list_of_triplets(genesDict[gene])[:-1] for gene in genesDict.keys()]
            codonList = [inner for outer in codonList for inner in outer]

            listy = cod_trna_aa_entropy.compute_compressibility_et_al([codonList], gc_model=True, prior_from_gc_model=True)
            metricWright.append(enc_wright(codonList))
            metricNovembre.append(enc_wright(codonList, novembreCorrection=True))
            metricInfoGC.append(listy[0]['overall_comp'])
            metricCodGivenTRNA.append(listy[1]['overall_comp'])
            metricTRNA.append(listy[2]['overall_comp'])
            genTimes.append(genomeDicty[x])

    print stats.spearmanr(metricWright, genTimes)
    print stats.spearmanr(metricNovembre, genTimes)
    print stats.spearmanr(metricInfoGC, genTimes)
    print stats.spearmanr(metricTRNA, genTimes)
    print stats.spearmanr(metricCodGivenTRNA, genTimes)


    return

def growth_metric_complicated():
    genomeDicty = {}
    f = open('../Data/growthRatesConvertedEdited.csv').readlines()
    for x in f:
        listy = x.split('\t')
        print listy
        listy[3] = listy[3].strip('\n')
        if len(listy[3]) > 5 and len(listy[2]) > 0:
            genomeDicty[listy[3].strip('\'').strip('"')] = float(listy[2])
    organismNames = []
    riboProteinGeneNumber = []
    allGenesNumber = []
    metricWrightAll = []
    metricNovembreAll = []
    metricSunAll = []
    metricInfoGCAll = []
    metricTRNAAll = []
    metricCodGivenTRNAAll = []
    metricWrightRibo = []
    metricNovembreRibo = []
    metricSunRibo = []
    metricInfoGCRibo = []
    metricTRNARibo = []
    metricCodGivenTRNARibo = []
    genTimes = []
    GCpercents = []
    outfile = open('../Results/GenerationTimeAnalysis.csv', 'w')
    for x in genomeDicty.keys()[:]:
        try:
            riboProtFile = '../Data/ribosomalProteins/' + x.split('/')[-2] + '.txt'
            f = open(riboProtFile).readlines()
            if len(f) > 5 and x.find('Mycoplasma') == -1:
                mycoplasma=False
            else:
                mycoplasma=True
            print '########################################################################' 
            print x
            print '########################################################################' 
            
            genesRibo = []
            lociNamesRibo = []
            for genbankFile in glob.glob(x + '*.gbk'):
                genesSingle, lociNamesSingle = pgb.get_coding_sequences(genbankFile, subsetFileLocation=riboProtFile, returnLociNames=True, mycoplasma=mycoplasma)
                if len(genesSingle) > len(genesRibo):
                    genesRibo = genesSingle
                    lociNamesRibo = lociNamesSingle
            genesDictRibo = {}
            for y in range(len(lociNamesRibo)):
                genesDictRibo[lociNamesRibo[y]] = genesRibo[y]
            print "I am using", len(genesDictRibo.keys()), "ribosomal genes for this organism"
            codonListRibo = [saf.split_string_into_list_of_triplets(genesDictRibo[gene])[:-1] for gene in genesDictRibo.keys()]
            codonListRibo = [inner for outer in codonListRibo for inner in outer]

            listyRibo = cod_trna_aa_entropy.compute_compressibility_et_al([codonListRibo], gc_model=True, prior_from_gc_model=True, mycoplasma=mycoplasma)
            metricWrightRibo.append(enc_wright(codonListRibo, mycoplasma=mycoplasma))
            metricNovembreRibo.append(enc_wright(codonListRibo, novembreCorrection=True, mycoplasma=mycoplasma))
            metricSunRibo.append(enc_sun(codonListRibo, mycoplasma=mycoplasma))
            metricInfoGCRibo.append(listyRibo[0]['overall_comp'])
            metricCodGivenTRNARibo.append(listyRibo[1]['overall_comp'])
            metricTRNARibo.append(listyRibo[2]['overall_comp'])
            
            genesAll = []
            lociNamesAll = []
            for genbankFile in glob.glob(x + '*.gbk'):
                genesSingle, lociNamesSingle = pgb.get_coding_sequences(genbankFile, returnLociNames=True, mycoplasma=mycoplasma)
                if len(genesSingle) > len(genesAll):
                    genesAll = genesSingle
                    lociNamesAll = lociNamesSingle
            genesDictAll = {}
            for y in range(len(lociNamesAll)):
                genesDictAll[lociNamesAll[y]] = genesAll[y]
            print "I am using", len(genesDictAll.keys()), "genes for this organism"
            
            codonListAll = [saf.split_string_into_list_of_triplets(genesDictAll[gene])[:-1] for gene in genesDictAll.keys()]
            codonListAll = [inner for outer in codonListAll for inner in outer]

            listyAll = cod_trna_aa_entropy.compute_compressibility_et_al([codonListAll], gc_model=True, prior_from_gc_model=True, mycoplasma=mycoplasma)
            metricWrightAll.append(enc_wright(codonListAll, mycoplasma=mycoplasma))
            metricNovembreAll.append(enc_wright(codonListAll, novembreCorrection=True, mycoplasma=mycoplasma))
            metricSunAll.append(enc_sun(codonListAll, mycoplasma=mycoplasma))
            metricInfoGCAll.append(listyAll[0]['overall_comp'])
            metricCodGivenTRNAAll.append(listyAll[1]['overall_comp'])
            metricTRNAAll.append(listyAll[2]['overall_comp'])
            riboProteinGeneNumber.append(len(genesDictRibo.keys())) 
            allGenesNumber.append(len(genesDictAll.keys()))
            organismNames.append(x)
            genTimes.append(genomeDicty[x])
            GCpercents.append(get_gc_percent(''.join(codonListAll)))
            outfile.write(x + ','\
                    + str(genomeDicty[x]) + ','\
                    + str(get_gc_percent(''.join(codonListAll))) + ','\
                    + str(len(genesDictRibo.keys())) + ','\
                    + str(len(genesDictAll.keys())) + ','\
                    + str(enc_wright(codonListRibo, mycoplasma=mycoplasma)) + ','\
                    + str(enc_wright(codonListAll, mycoplasma=mycoplasma)) + ','\
                    + str(enc_wright(codonListRibo, novembreCorrection=True, mycoplasma=mycoplasma)) + ','\
                    + str(enc_wright(codonListAll, novembreCorrection=True, mycoplasma=mycoplasma)) + ','\
                    + str(enc_sun(codonListRibo, mycoplasma=mycoplasma)) + ','\
                    + str(enc_sun(codonListAll, mycoplasma=mycoplasma)) + ','\
                    + str(listyRibo[0]['overall_comp']) + ','\
                    + str(listyAll[0]['overall_comp'])\
                    + '\n')

        except IOError:
            pass
    outfile.close()
    

    
    metricWright = [(metricWrightAll[x] - metricWrightRibo[x])/float(metricWrightAll[x]) for x in range(len(genTimes))]
    metricNovembre = [(metricNovembreAll[x] - metricNovembreRibo[x])/float(metricNovembreAll[x]) for x in range(len(genTimes))]
    metricSun = [(metricSunAll[x] - metricSunRibo[x])/float(metricSunAll[x]) for x in range(len(genTimes))]
    metricInfoGC = [(metricInfoGCAll[x] - metricInfoGCRibo[x])/float(metricInfoGCAll[x]) for x in range(len(genTimes))]
    metricTRNA = [(metricTRNAAll[x] - metricTRNARibo[x])/float(metricTRNAAll[x]) for x in range(len(genTimes))]
    metricCodGivenTRNA = [(metricCodGivenTRNAAll[x] - metricCodGivenTRNARibo[x])/float(metricCodGivenTRNAAll[x]) for x in range(len(genTimes))]
    
    
    
    print "The total Analysis included", len(genTimes), "genomes"
    print "Wright"
    print stats.spearmanr(metricWright, genTimes)
    print stats.spearmanr(metricWrightAll, genTimes)
    print stats.spearmanr(metricWrightRibo, genTimes)
    print "Novembre (published)"
    print stats.spearmanr(metricNovembre, genTimes)
    print stats.spearmanr(metricNovembreAll, genTimes)
    print stats.spearmanr(metricNovembreRibo, genTimes)
    print "Sun"
    print stats.spearmanr(metricSun, genTimes)
    print stats.spearmanr(metricSunAll, genTimes)
    print stats.spearmanr(metricSunRibo, genTimes)
    print "InfoBias: Codon given AA"
    print stats.spearmanr(metricInfoGC, genTimes)
    print stats.spearmanr(metricInfoGCAll, genTimes)
    print stats.spearmanr(metricInfoGCRibo, genTimes)
    print "InfoBias: tRNA given AA"
    print stats.spearmanr(metricTRNA, genTimes)
    print stats.spearmanr(metricTRNAAll, genTimes)
    print stats.spearmanr(metricTRNARibo, genTimes)
    print "InfoBias: Codon given tRNA"
    print stats.spearmanr(metricCodGivenTRNA, genTimes)
    print stats.spearmanr(metricCodGivenTRNAAll, genTimes)
    print stats.spearmanr(metricCodGivenTRNARibo, genTimes)

    print "Correlations with GC percent"
    print stats.spearmanr(GCpercents, genTimes)
    print "Wright"
    print stats.spearmanr(metricWright, GCpercents)
    print stats.spearmanr(metricWrightAll, GCpercents)
    print stats.spearmanr(metricWrightRibo, GCpercents)
    print "Novembre (published)"
    print stats.spearmanr(metricNovembre, GCpercents)
    print stats.spearmanr(metricNovembreAll, GCpercents)
    print stats.spearmanr(metricNovembreRibo, GCpercents)
    print "Sun"
    print stats.spearmanr(metricSun, GCpercents)
    print stats.spearmanr(metricSunAll, GCpercents)
    print stats.spearmanr(metricSunRibo, GCpercents)
    print "InfoBias: Codon given AA"
    print stats.spearmanr(metricInfoGC, GCpercents)
    print stats.spearmanr(metricInfoGCAll, GCpercents)
    print stats.spearmanr(metricInfoGCRibo, GCpercents)
    


    return

def enc_wright_aa(tempCountsDict, aaDict, aminoAcid, prior=False):
    #n is the total occurrences of the amino acid
    n = float(sum(tempCountsDict.values()))
    k = float(len(aaDict[aminoAcid]))
    #need to have n > 1 for the amino acid for anything to make any sense
    if n > 5:
        tempProbs = {}
        for codon in aaDict[aminoAcid]:
            tempProbs[codon] = ((tempCountsDict[codon]/n)**2)
        F = ((n * sum(tempProbs.values()))-1)/(n-1)
        #Put all F values into the dictionary for use later
    #if n<=1 we're going to throw it away
    elif n <= 5:
        F = 'n/a'
    return F

def enc_novembre_aa(tempCountsDict, aaDict, aminoAcid, nucFreqsDict):
    #n is the total occurrences of the amino acid
    n = float(sum(tempCountsDict.values()))
    #need to have n > 1 for the amino acid for anything to make any sense
    k = len(aaDict[aminoAcid])
    if n > 5:
        expectedDict = {}
        for codon in aaDict[aminoAcid]:
            first = codon[0]
            second = codon[1]
            third = codon[2]
            expectedDict[codon] = nucFreqsDict[first]*nucFreqsDict[second]*nucFreqsDict[third]
        reNormConstant = 1./sum(expectedDict.values())
        for codon in aaDict[aminoAcid]:
            expectedDict[codon] = expectedDict[codon] * reNormConstant
        tempProbs = {}
        for codon in aaDict[aminoAcid]:
            tempProbs[codon] = ((tempCountsDict[codon]/n))
            tempProbs[codon] = (n*((tempProbs[codon] - expectedDict[codon])**2))/float(expectedDict[codon])
        chiSquared = sum(tempProbs.values())
        F = (chiSquared+n-k) / (k*(n-1))
        #Put all F values into the dictionary for use later
    #if n<=1 we're going to throw it away
    elif n <= 5:
        F = 'n/a'
    return F

def enc_sun_aa(tempCountsDict, aaDict, aminoAcid, prior):
    #n is the total occurrences of the amino acid
    n = float(sum(tempCountsDict.values()))
    tempProbs = {}
    for codon in aaDict[aminoAcid]:
        #essentially this is a squared probability while adding in pseudo counts
        try:
            tempProbs[codon] = ((tempCountsDict[codon]+prior)/(n+(prior*len(aaDict[aminoAcid]))))
        except ZeroDivisionError:
            tempProbs[codon] = 0

    Fcf = sum([x**2 for x in tempProbs.values()])
    n = n + (prior*len(aaDict[aminoAcid]))
    return Fcf, n


########################################################################################
#The classical and new version of the Effective Number of Codons
def enc_wright(codonList, novembreCorrection=False, mycoplasma=False):
    '''Takes a list of codons and outputs the value of the Effective Number of Codons
    as calculated by Wright 1990. Uses a lot of hacks because the original formulation 
    itself is quite hacky'''
    if mycoplasma:
        aaDict = saf.get_AA_dict(mycoplasma=True)
    else:
        aaDict = saf.get_AA_dict()
    aaEncDict = {}
    
    if novembreCorrection:
        nucFreqsDict = {}
        codonString = ''.join(codonList)
        for nuc in ['A', 'T', 'G', 'C']:
            nucFreqsDict[nuc] = codonString.count(nuc)/float(len(codonString))

    #iterate through each amino  acid
    for aminoAcid in aaDict.keys():
        tempCounts = {}
        for codon in aaDict[aminoAcid]:
            #iterate through each codon for that amino acid and count occurrences in codonList
            tempCounts[codon] = codonList.count(codon)
        if novembreCorrection:
            F = enc_novembre_aa(tempCounts, aaDict, aminoAcid, nucFreqsDict)
        else:
            F = enc_wright_aa(tempCounts, aaDict, aminoAcid)
        aaEncDict[aminoAcid] = F

    #There are 4 classes of amino acids depending on how many codons they have that code for each
    #They will be treated separately here.
    #First separate the F values accordingly
    twos = [aaEncDict[aminoAcid] for aminoAcid in aaDict.keys() if len(aaDict[aminoAcid]) == 2]
    #And remove any F values that were a string, i.e. n/a
    twos = [x for x in twos if not isinstance(x, str)]
    #Repeat process for three-fold redundant codon families
    threes = [aaEncDict[aminoAcid] for aminoAcid in aaDict.keys() if len(aaDict[aminoAcid]) == 3]
    threes = [x for x in threes if not isinstance(x, str)]
    #And fours
    fours = [aaEncDict[aminoAcid] for aminoAcid in aaDict.keys() if len(aaDict[aminoAcid]) == 4]
    fours = [x for x in fours if not isinstance(x, str)]
    #And sixes
    sixes = [aaEncDict[aminoAcid] for aminoAcid in aaDict.keys() if len(aaDict[aminoAcid]) == 6]
    sixes = [x for x in sixes if not isinstance(x, str)]
    
    ##########################################################################
    #Test whether we can actually do anything here and return nothing if we can't
    if len(twos) == 0 or np.mean(twos) < 1/float(2):
        twos = [1/float(2)]
    elif len(fours) == 0 or np.mean(fours) < 1/float(4):
        fours = [1/float(4)]
    elif len(sixes) == 0 or np.mean(fours) < 1/float(6):
        sixes = [1/float(6)]

    #Onto the actual calculation. 
    twosF = 9./float(np.mean(twos))
    foursF = 5./float(np.mean(fours))
    sixesF = 3./float(np.mean(sixes))
    #There is only one threefold redundant family so it will be treated differently
    #Namely, if there is not enough of it to actually calculate we will use an average
    #from the twofold and fourfold redundant families to approximate its value
    if len(threes) == 1:
        try:
            threesF = 1./float(threes[0])
        except ZeroDivisionError:
            threesF = 1./float(np.mean([np.mean(twos), np.mean(fours)])) 
    elif len(threes) == 0:
        threesF = 1./float(np.mean([np.mean(twos), np.mean(fours)])) 
    elif len(threes) > 1:
        print "Big problem here, why is there more than one amino acid with three codons?"
        return 'n/a'
    ENC = 2 + twosF + threesF + foursF + sixesF
#how to deal with the fact that ENC's often get over 61? a hacky way is to just throw them all back to 61
    return ENC

def enc_sun(codonList, prior=1, mycoplasma=False):
    '''
    Takes a list of codons and outputs the value of the Effective Number of Codons
    as updated by Sun, Yang and Xia 2012.
    One difference is that I am not separating the 6-fold amino acids into 2 groups. Because why would I?
    '''
    
    if mycoplasma:
        aaDict = saf.get_AA_dict(mycoplasma=True)
    else:
        aaDict = saf.get_AA_dict()
    aaDict['leucine'] = ['CTT', 'CTC', 'CTA', 'CTG']
    aaDict['arginine'] = ['CGT', 'CGC', 'CGA', 'CGG']
    aaDict['serine'] = ['TCT', 'TCC', 'TCA', 'TCG']
    aaDict['leucine2'] = ['TTA', 'TTG']
    aaDict['arginine2'] = ['AGA', 'AGG']
    aaDict['serine2'] = ['AGT', 'AGC']

    aaEncDict = {}
    aaNDict = {}
    #iterate through each amino acid
    for aminoAcid in aaDict.keys():
        tempCounts = {}
        for codon in aaDict[aminoAcid]:
            #iterate through each codon for that amino acid and count occurrences in codonList
            tempCounts[codon] = codonList.count(codon)
        aaEncDict[aminoAcid], aaNDict[aminoAcid] = enc_sun_aa(tempCounts, aaDict, aminoAcid, prior=prior) 
    #There are 3 classes of amino acids depending on how many codons they have that code for each
    #They will be treated separately here.
    #First separate the F and N values accordingly
    twosF = [aaEncDict[aminoAcid] for aminoAcid in aaDict.keys() if len(aaDict[aminoAcid]) == 2]
    twosN = [aaNDict[aminoAcid] for aminoAcid in aaDict.keys() if len(aaDict[aminoAcid]) == 2]
    #Repeat process for three-fold redundant codon families
    threesF = [aaEncDict[aminoAcid] for aminoAcid in aaDict.keys() if len(aaDict[aminoAcid]) == 3]
    threesN = [aaNDict[aminoAcid] for aminoAcid in aaDict.keys() if len(aaDict[aminoAcid]) == 3]
    #And fours
    foursF = [aaEncDict[aminoAcid] for aminoAcid in aaDict.keys() if len(aaDict[aminoAcid]) == 4]
    foursN = [aaNDict[aminoAcid] for aminoAcid in aaDict.keys() if len(aaDict[aminoAcid]) == 4]
    #And sixes
    sixesF = [aaEncDict[aminoAcid] for aminoAcid in aaDict.keys() if len(aaDict[aminoAcid]) == 6]
    sixesN = [aaNDict[aminoAcid] for aminoAcid in aaDict.keys() if len(aaDict[aminoAcid]) == 6]
    ############################################################################################33
    #now take the weighted average
    try:
        twos = (len(twosN) * sum(twosN))/float(sum([x*y for x,y in zip(twosF, twosN)]))
    except ZeroDivisionError:
        twos = 0
    try:
        threes = (len(threesF) * sum(threesN))/float(sum([x*y for x,y in zip(threesF, threesN)]))
    except ZeroDivisionError:
        threes = 0
    try:
        fours = (len(foursF) * sum(foursN))/float(sum([x*y for x,y in zip(foursF, foursN)]))
    except ZeroDivisionError:
        fours = 0
    try:
        sixes = (len(sixesF) * sum(sixesN))/float(sum([x*y for x,y in zip(sixesF, sixesN)]))
    except ZeroDivisionError:
        sixes = 0
    ENC = 2 + twos + threes + fours + sixes
    return ENC

def get_gc_percent(geneString):
    GC = 0
    for x in geneString:
        if x in ['G', 'C']:
            GC+=1
    GC = GC/float(len(geneString))
    return GC

########################################################################################
########################################################################################
#Plots
########################################################################################
########################################################################################
def read_paxdb(expressionFile):
    dicty = {}
    inFile = open(expressionFile, 'r').readlines()
    for line in inFile:
        if line[0] != '#':
            splitLine = line.split('\t')
            dicty[splitLine[1].split('.')[1]] = float(splitLine[2])
    return dicty

def all_other_organism_expression():
#    organisms = ['escherichia', 'leptospira', 'mycobacterium', 'mycoplasma', 'shigellaF', 'shigellaD', 'bacillus', 'halobacterium', 'thermococcus', 'salmonella', 'streptococcus', 'microcystis', 'synechocystis', 'deinococcus']
    organisms = ['escherichia', 'leptospira', 'mycobacterium', 'shigellaF', 'shigellaD', 'bacillus', 'halobacterium', 'thermococcus', 'salmonella', 'streptococcus', 'microcystis', 'synechocystis', 'deinococcus']
    dicty = {}
    problemOrganisms = []
    for organism in organisms:
        if organism == 'mycoplasma':
            a, b, c, d, e = other_organism_expression(organism, mycoplasma=True)
        else: 
            a, b, c, d, e = other_organism_expression(organism)

        dicty[organism] = [a,b,c,d,e]
    return dicty

def other_organism_expression_plot(dicty):    
    Original = []
    Novembre = []
    Sun = []
    GCprior = []
    for organism in dicty.keys():
        Original.append(dicty[organism][0][0])
        Novembre.append(dicty[organism][1][0])
        Sun.append(dicty[organism][2][0])
        GCprior.append(dicty[organism][4][0])
    
    print stats.wilcoxon(Original, GCprior)
    print stats.wilcoxon(Novembre, GCprior)
    print stats.wilcoxon(Sun, GCprior)
    print stats.wilcoxon(Original, Novembre)
    
    print stats.ttest_rel(Original, GCprior)
    print stats.ttest_rel(Novembre, GCprior)
    print stats.ttest_rel(Sun, GCprior)
    print stats.ttest_rel(Original, Novembre)

    #plt.close()
    #fig = plt.figure()
    #plt.boxplot([Original, Novembre, Sun, GCprior])
    #plt.tight_layout()
    #repeats = 6
    #colors = ['b', 'g', 'r', 'c', 'y', 'k']
    #colors=colors*repeats
    #tempAlphas = [1, 0.9, 0.8, 0.7, 0.6, 0.5] 
    #alphas = [tempAlphas[y] for y in range(repeats) for x in range(7)]
    #beeswarm.beeswarm([Original, Novembre, Sun, GCprior], col=colors[:len(Original)], method='center', s=400) 
    #plt.show()

    repeats = 4
    xData = [0,1,2,3]
    offsetX = [0.05, 0.95]
    pairs = [(p,q,r,s) for (p,q,r,s) in zip(Original,Novembre,Sun,GCprior)]
    colors = ['b', 'g', 'r', 'c', 'm', 'y', 'k']
    colors = colors*repeats
    tempAlphas = [0.9, 0.5, 0.4, 0.3]
    alphas = [tempAlphas[y] for y in range(repeats) for x in range(7)]
    for data in range(len(pairs)):
        print colors[data], alphas[data]
        plt.scatter(xData, pairs[data], c=colors[data], alpha=alphas[data], edgecolor='none', s=200)
        #plt.plot(offsetX, pairs[data], colors[data]+'--', alpha=alphas[data])
    plt.show()
    return

def other_organism_expression(speciesName, mycoplasma=False):
    if speciesName == 'escherichia':
        genbankFiles = ['../../GenBank/allBacteria/Escherichia_coli_K_12_substr__MG1655_uid57779/NC_000913.gbk']
        expressionFile = '../Data/expressionData/escherichiaIntegrated.txt'
    
    elif speciesName == 'leptospira':
        genbankFiles = ['../../GenBank/allBacteria/Leptospira_interrogans_serovar_Copenhageni_Fiocruz_L1_130_uid58065/NC_005823.gbk',\
                '../../GenBank/allBacteria/Leptospira_interrogans_serovar_Copenhageni_Fiocruz_L1_130_uid58065/NC_005824.gbk']
        expressionFile = '../Data/expressionData/leptospiraIntegrated.txt'
    
    elif speciesName == 'mycobacterium':
        genbankFiles = ['../../GenBank/allBacteria/Mycobacterium_tuberculosis_H37Rv_uid57777/NC_000962.gbk']
        expressionFile = '../Data/expressionData/mycobacteriumIntegrated.txt'
    
    elif speciesName == 'mycoplasma':
        genbankFiles = ['../../GenBank/allBacteria/Mycoplasma_pneumoniae_M129_uid57709/NC_000912.gbk']
        expressionFile = '../Data/expressionData/mycoplasma.txt'
    
    elif speciesName == 'shigellaF':
        genbankFiles = ['../../GenBank/allBacteria/Shigella_flexneri_2a_2457T_uid57991/NC_004741.gbk']
        expressionFile = '../Data/expressionData/shigellaF.txt'
    
    elif speciesName == 'shigellaD':
        genbankFiles = ['../../GenBank/allBacteria/Shigella_dysenteriae_Sd197_uid58213/NC_007606.gbk',\
                '../../GenBank/allBacteria/Shigella_dysenteriae_Sd197_uid58213/NC_009344.gbk',\
                '../../GenBank/allBacteria/Shigella_dysenteriae_Sd197_uid58213/NC_007607.gbk']
        expressionFile = '../Data/expressionData/shigellaD.txt'
    
    elif speciesName == 'bacillus':
        genbankFiles = ['../../GenBank/allBacteria/Bacillus_subtilis_168_uid57675/NC_000964.gbk']
        expressionFile = '../Data/expressionData/bacillus.txt'
    
    elif speciesName == 'halobacterium':
        genbankFiles = ['../../GenBank/allBacteria/Halobacterium_NRC_1_uid57769/NC_001869.gbk',\
                '../../GenBank/allBacteria/Halobacterium_NRC_1_uid57769/NC_002607.gbk',\
                '../../GenBank/allBacteria/Halobacterium_NRC_1_uid57769/NC_002608.gbk']
        expressionFile = '../Data/expressionData/halobacterium.txt'
    
    elif speciesName == 'thermococcus':
        genbankFiles = ['../../GenBank/allBacteria/Thermococcus_gammatolerans_EJ3_uid59389/NC_012804.gbk']
        expressionFile = '../Data/expressionData/thermococcus.txt'
    
    elif speciesName == 'salmonella':
        genbankFiles = ['../../GenBank/allBacteria/Salmonella_enterica_serovar_Typhimurium_LT2_uid57799/NC_003197.gbk',\
                '../../GenBank/allBacteria/Salmonella_enterica_serovar_Typhimurium_LT2_uid57799/NC_003277.gbk']
        expressionFile = '../Data/expressionData/salmonella.txt'
    
    elif speciesName == 'streptococcus':
        genbankFiles = ['../../GenBank/allBacteria/Streptococcus_pyogenes_M1_GAS_uid57845/NC_002737.gbk']
        expressionFile = '../Data/expressionData/streptococcusIntegrated.txt'
    
    elif speciesName == 'microcystis':
        genbankFiles = ['../../GenBank/allBacteria/Microcystis_aeruginosa_NIES_843_uid59101/NC_010296.gbk']
        expressionFile = '../Data/expressionData/microcystisIntegrated.txt'
    
    elif speciesName == 'synechocystis':
        genbankFiles = ['../../GenBank/allBacteria/Synechocystis_PCC_6803_uid57659/NC_000911.gbk',\
                '../../GenBank/allBacteria/Synechocystis_PCC_6803_uid57659/NC_005229.gbk',\
                '../../GenBank/allBacteria/Synechocystis_PCC_6803_uid57659/NC_005230.gbk',\
                '../../GenBank/allBacteria/Synechocystis_PCC_6803_uid57659/NC_005231.gbk',\
                '../../GenBank/allBacteria/Synechocystis_PCC_6803_uid57659/NC_005232.gbk']
        expressionFile = '../Data/expressionData/synechocystis.txt'
    
    elif speciesName == 'deinococcus':
        genbankFiles = ['../../GenBank/allBacteria/Deinococcus_deserti_VCD115_uid58615/NC_012526.gbk',\
                '../../GenBank/allBacteria/Deinococcus_deserti_VCD115_uid58615/NC_012527.gbk',\
                '../../GenBank/allBacteria/Deinococcus_deserti_VCD115_uid58615/NC_012528.gbk',\
                '../../GenBank/allBacteria/Deinococcus_deserti_VCD115_uid58615/NC_012529.gbk']
        expressionFile = '../Data/expressionData/deinococcus.txt'


    genes = []
    lociNames = []
    for genbankFile in range(len(genbankFiles)):
        if speciesName != 'mycoplasma':
            genesSingle, lociNamesSingle = pgb.get_coding_sequences(genbankFiles[genbankFile], returnLociNames=True)
        else:
            genesSingle, lociNamesSingle = pgb.get_coding_sequences(genbankFiles[genbankFile], returnLociNames=True, mycoplasma=True)
        for i in range(len(genesSingle)):
            genes.append(genesSingle[i])
            lociNames.append(lociNamesSingle[i])
    
    #Put them into a dictionary
    genesDict = {}
    for x in range(len(lociNames)):
        genesDict[lociNames[x]] = genes[x]

    proteinDictionary = read_paxdb(expressionFile)
    
    GCcontent = get_gc_percent(''.join(genes))

    encPrior = []
    encGCPrior = []
    length = []
    protein = []
    wrightENC = []
    novembreENC = []
    sunENC = []
    
    for lociName in genesDict.keys():
        if lociName in proteinDictionary.keys():
            codonList = saf.split_string_into_list_of_triplets(genesDict[lociName])[:-1]
            #listy = cod_trna_aa_entropy.compute_compressibility_et_al([codonList], gc_model=False, prior_aa=1, mycoplasma=mycoplasma)
            #encPrior.append(listy[0]['overall_comp'])
            #listy = cod_trna_aa_entropy.compute_compressibility_et_al([codonList], gc_model=True, prior_from_gc_model=True, mycoplasma=mycoplasma)
            #encGCPrior.append(listy[0]['overall_comp'])
            listy = cod_trna_aa_entropy.compute_compressibility_et_al([codonList], gc_model=False, prior_aa=1)
            encPrior.append(listy[1]['overall_comp'])
            listy = cod_trna_aa_entropy.compute_compressibility_et_al([codonList], gc_model=True, prior_from_gc_model=True)
            encGCPrior.append(listy[1]['overall_comp'])
            
            length.append(len(codonList))
            protein.append(proteinDictionary[lociName])
            wrightENC.append(enc_wright(codonList, mycoplasma=mycoplasma))
            novembreENC.append(enc_wright(codonList, novembreCorrection=True, mycoplasma=mycoplasma))
            sunENC.append(enc_sun(codonList, prior=1, mycoplasma=mycoplasma))
    
    print "Species:", speciesName, "has", len(length), "data points and a GC content of", GCcontent
    print "The correlation between expression and gene length is:"
    lCorr = stats.spearmanr(length, protein)
    print lCorr

    print "Our non-GC corrected model with a prior a.a. of 1:"
    noGCyesPriorCorr = stats.spearmanr(protein, encPrior)
    print noGCyesPriorCorr
    
    print "Our GC corrected model with a prior a.a. of 1:"
    yesGCyesPriorCorr = stats.spearmanr(protein, encGCPrior)
    print yesGCyesPriorCorr

    print "For comparison, the original formulation correlation coefficient is:"
    originalCorr = stats.spearmanr(protein, wrightENC)
    print originalCorr
    
    print "And the GC corrected version from Novembre is:"
    novembreCorr = stats.spearmanr(protein, novembreENC)
    print novembreCorr
    
    print "Whereas the methodology of Sun et al produces:"
    sunCorr = stats.spearmanr(protein, sunENC)
    print sunCorr

    return originalCorr, novembreCorr, sunCorr, noGCyesPriorCorr, yesGCyesPriorCorr

def enc_and_compress_for_ind_aas(specifiedSequence=False):
    '''
    By default, this will show a bar graph of the Effective number of codons for each
    individual amino acid for the genome specified below
    '''
    if not specifiedSequence:
        #Just getting the sequences to use from either a whole genome or a reference set 
        genbankFile = '../../CodonAnalysis/data/Genomes/escherichia/escherichia.gb'
#        genbankFile = '../../CodonAnalysis/data/Genomes/helicobacter/helicobacter.gb'
#        genbankFile = '../../CodonAnalysis/data/Genomes/pseudomonas/pseudomonas.gb'
#        genbankFile = '../../CodonAnalysis/data/Genomes/bacillus/bacillus.gb'
#        refSetFile = '../../CodonAnalysis/data/referenceSets/escherichiaBottom50percent'
#        refSetFile = '../../CodonAnalysis/data/referenceSets/escherichiaOriginalCAIset'
#        genes = pgb.get_coding_sequences(genbankFile, subsetFileLocation=refSetFile)
        genes, lociNames = pgb.get_coding_sequences(genbankFile, returnLociNames=True)
        #Put them into a dictionary
        dicty = {}
        for x in range(len(lociNames)):
            dicty[lociNames[x]] = genes[x]
    else:
        genes = [specifiedSequence]

    #Load a few basic dictionaries to deal with the data
    aaDict = saf.get_AA_dict(stop=False)
    codonToAAdict = saf.get_codon_to_AA_dict(stop=False)
    allCodons = aaDict.values()
    allCodons = [inner for outer in allCodons for inner in outer] 
    
    ###Creates the type of data structure I need which is a giant list of codons
    codonList = [saf.split_string_into_list_of_triplets(gene)[:-1] for gene in genes]
    codonList = [inner for outer in codonList for inner in outer]
    #################################################################################
    #main function
    #################################################################################
    listy = cod_trna_aa_entropy.compute_compressibility_et_al([codonList], gc_model=False)
    aaCodonDict = listy[0]
    trnaCodonDict = listy[1]
    aaTRNAdict = listy[2]

    #Get simple lists of the amino acids that have x number of codons
    sixList = []
    fourList = []
    threeList = []
    twoList = []
    oneList = []
    for aa in aaDict.keys():
        if len(aaDict[aa]) == 6:
            sixList.append(aa)
        if len(aaDict[aa]) == 4:
            fourList.append(aa)
        if len(aaDict[aa]) == 3:
            threeList.append(aa)
        if len(aaDict[aa]) == 2:
            twoList.append(aa)
        if len(aaDict[aa]) == 1:
            oneList.append(aa)
    
    #Split my data up according to the above groupings
    #NOTE: this section dictates which data I am caring about. i.e. codon given aa, trna given aa, etc.
    #as well as whether to care about Enc (1) or compressibility (0)
    geneAAdict6 = {}
    geneAAdict4 = {}
    geneAAdict3 = {}
    geneAAdict2 = {}
    geneAAdict1 = {}
    for x in aaCodonDict['x_to_comp_effno'].keys():
        if x in sixList:
            geneAAdict6[x] = aaCodonDict['x_to_comp_effno'][x][0]
        if x in fourList:
            geneAAdict4[x] = aaCodonDict['x_to_comp_effno'][x][0]
        if x in threeList:
            geneAAdict3[x] = aaCodonDict['x_to_comp_effno'][x][0]
        if x in twoList:
            geneAAdict2[x] = aaCodonDict['x_to_comp_effno'][x][0]
        if x in oneList:
            geneAAdict1[x] = aaCodonDict['x_to_comp_effno'][x][0]
    six = len(geneAAdict6)
    four = len(geneAAdict4)
    three = len(geneAAdict3)
    two = len(geneAAdict2)
    one = len(geneAAdict1)
   
    plt.figure() 
    plt.axvline(1, color='k')
    plt.barh(range(one,one+two), geneAAdict2.values(), align='center', facecolor='k', alpha=0.25)
    plt.barh(range(one+two, one+two+three), geneAAdict3.values(), align='center', facecolor='k', alpha=0.75)
    plt.barh(range(one+two+three, one+two+three+four), geneAAdict4.values(), align='center', facecolor='k', alpha=0.5)
    plt.barh(range(one+two+three+four, one+two+three+four+six), geneAAdict6.values(), align='center', facecolor='k', alpha=1)
    plt.yticks(range(one+two+three+four+six), geneAAdict1.keys()+geneAAdict2.keys()+geneAAdict3.keys()+geneAAdict4.keys()+geneAAdict6.keys(), ha='center')
    plt.xlim(0,1.1)
    plt.ylim(-1,18)
    plt.tick_params(labelsize=20)
    plt.tick_params(axis='y', pad=60)
    plt.tight_layout()

    geneAAdict6 = {}
    geneAAdict4 = {}
    geneAAdict3 = {}
    geneAAdict2 = {}
    geneAAdict1 = {}
    for x in aaCodonDict['x_to_comp_effno'].keys():
        if x in sixList:
            geneAAdict6[x] = aaCodonDict['x_to_comp_effno'][x][1]
        if x in fourList:
            geneAAdict4[x] = aaCodonDict['x_to_comp_effno'][x][1]
        if x in threeList:
            geneAAdict3[x] = aaCodonDict['x_to_comp_effno'][x][1]
        if x in twoList:
            geneAAdict2[x] = aaCodonDict['x_to_comp_effno'][x][1]
        if x in oneList:
            geneAAdict1[x] = aaCodonDict['x_to_comp_effno'][x][1]
    six = len(geneAAdict6)
    four = len(geneAAdict4)
    three = len(geneAAdict3)
    two = len(geneAAdict2)
    one = len(geneAAdict1)
    plt.figure() 
    plt.axvline(2, 0, 0.5, color='k')
    plt.axvline(3, 0, 0.55, color='k')
    plt.axvline(4, 0, 0.81, color='k')
    plt.axvline(6, 0, 1, color='k')
    plt.barh(range(one,one+two), geneAAdict2.values(), align='center', facecolor='k', alpha=0.25)
    plt.barh(range(one+two, one+two+three), geneAAdict3.values(), align='center', facecolor='k', alpha=0.75)
    plt.barh(range(one+two+three, one+two+three+four), geneAAdict4.values(), align='center', facecolor='k', alpha=0.5)
    plt.barh(range(one+two+three+four, one+two+three+four+six), geneAAdict6.values(), align='center', facecolor='k', alpha=1)
    #plt.yticks(range(one+two+three+four+six), geneAAdict1.keys()+geneAAdict2.keys()+geneAAdict3.keys()+geneAAdict4.keys()+geneAAdict6.keys())
    plt.yticks(range(one+two+three+four+six), [''])
    plt.xlim(0,6.2)
    plt.ylim(-1,18)
    plt.tick_params(labelsize=20)
    plt.tight_layout()

    
    plt.show()
    return 

def amino_acid_heatmaps(sorting='PNAS-digitalCountsAvg'):
    '''
    Heatmaps of compressibility for different models
    '''
    #Just getting the sequences to use from either a whole genome or a reference set 
    genbankFile = '../../CodonAnalysis/data/Genomes/escherichia/escherichia.gb'
#    genbankFile = '../../CodonAnalysis/data/Genomes/helicobacter/helicobacter.gb'
#    refSetFile = '../../CodonAnalysis/data/referenceSets/escherichiaBottom50percent'
#    refSetFile = '../../CodonAnalysis/data/referenceSets/escherichiaOriginalCAIset'
#    genes = pgb.get_coding_sequences(genbankFile, subsetFileLocation=refSetFile)
    genes, lociNames = pgb.get_coding_sequences(genbankFile, returnLociNames=True)
    #Put them into a dictionary
    genesDict = {}
    for x in range(len(lociNames)):
        genesDict[lociNames[x]] = genes[x]

    ###############################################
    #loading some basic dictionaries
    aaDict = saf.get_AA_dict(stop=False)
    aaDictNew = saf.get_AA_dict(stop=False)
    for aa in aaDictNew.keys():
        aaDictNew[aa] = 0
    codonToAAdict = saf.get_codon_to_AA_dict(stop=False)
    allCodons = codonToAAdict.keys()
    trnaDict = get_trna_dict()
    
    allGenesCodonAA = {}
    allGenesCodonTRNA = {}
    allGenesTRNAaa = {}
    allGenesJustAA = {}
    allGenesJustCodon = {}
    overallCodonGivenAA = {}
    overallTRNAgivenAA = {}
    overallCodonGivenTRNA = {}
    for lociName in genesDict.keys():
        codonList = saf.split_string_into_list_of_triplets(genesDict[lociName])[:-1]
        listy = cod_trna_aa_entropy.compute_compressibility_et_al([codonList], gc_model=True, prior_from_gc_model=True)
#        listy = cod_trna_aa_entropy.compute_compressibility_et_al([codonList], gc_model=False, prior_aa=1)
        aaCodonDict = listy[0]
        trnaCodonDict = listy[1]
        aaTRNAdict = listy[2]

        allGenesCodonAA[lociName] = split_tuple_dict(aaDict, aaCodonDict['x_to_comp_effno'])
        allGenesCodonTRNA[lociName] = split_tuple_dict(trnaDict, trnaCodonDict['x_to_comp_effno'])
        allGenesTRNAaa[lociName] = split_tuple_dict(aaDict, aaTRNAdict['x_to_comp_effno'])
        allGenesJustAA[lociName] = aaCodonDict['eff_x']
        allGenesJustCodon[lociName] = aaCodonDict['eff_y']
        
        overallTRNAgivenAA[lociName] = aaTRNAdict['overall_comp']
        overallCodonGivenTRNA[lociName] = trnaCodonDict['overall_comp']
        overallCodonGivenAA[lociName] = aaCodonDict['overall_comp']
    
    ####using expression to sort
    expressionDict = compileExpressionData.read_protein_expression_data(sorting, '../../GeneralFunctions/Data/proteinExpressionP.csv')
    toDelete = []
    for key in expressionDict.keys():
        if key not in genesDict.keys():
            del expressionDict[key]

    sortedData = sorted(expressionDict.iteritems(), key=operator.itemgetter(1), reverse=True)
    sortedLociNames = [x[0] for x in sortedData]

    listA = [overallTRNAgivenAA[lociName] for lociName in expressionDict.keys()]
    listB = [overallCodonGivenTRNA[lociName] for lociName in expressionDict.keys()]
    listC = [overallCodonGivenAA[lociName] for lociName in expressionDict.keys()]
    expressionValues = [expressionDict[x][0] for x in expressionDict.keys()]
    
    print "Correlation with expression data:"
    print "codon given aa", stats.spearmanr(listC, expressionValues)
    print "codon given trna", stats.spearmanr(listB, expressionValues)
    print "trna given a", stats.spearmanr(listA, expressionValues)
    

    ###################################################################################
    #plotting the data 
    ###################################################################################
    plt.close()

    ###################################################################################
    #the effective number of amino acids
    allGenesList = []
    for x in sortedLociNames:
        try:
            allGenesList.append([allGenesJustAA[x]])
        except KeyError:
            pass
    
    fig = plt.figure()
    ax1 = fig.add_subplot(111)
    results = ax1.matshow(np.array(allGenesList), extent = [0, len(allGenesList[0]), 0, len(allGenesList)], aspect='auto')
#    plt.xticks(range(0, len(aminoAcidOrder), 1), rotation=90)
#    ax1.set_xticklabels(aminoAcidOrder)
    plt.colorbar(results)


    ###################################################################################     
    #codon given AA
    aminoAcidOrder = get_aa_sort_list(genesDict, typeFlag='aaCodon')
    allGenesList = []
    for x in sortedLociNames:
        try:
            allGenesList.append(sort_amino_acids(allGenesCodonAA[x], aminoAcidOrder))
        except KeyError:
            pass
    fig = plt.figure()
    ax1 = fig.add_subplot(111)
    results = ax1.matshow(np.array(allGenesList), extent = [0, len(allGenesList[0]), 0, len(allGenesList)], aspect='auto')
    plt.xticks(range(0, len(aminoAcidOrder), 1), rotation=90)
    ax1.set_xticklabels(aminoAcidOrder, ha='left')
    plt.colorbar(results)
    plt.tight_layout()

    ###################################################################################     
    #codon given trna 
    aminoAcidOrder = get_aa_sort_list(genesDict, typeFlag='trnaCodon')
    allGenesList = []
    for x in sortedLociNames:
        try:
            allGenesList.append(sort_amino_acids(allGenesCodonTRNA[x], aminoAcidOrder))
        except KeyError:
            pass
    fig = plt.figure()
    ax1 = fig.add_subplot(111)
    results = ax1.matshow(np.array(allGenesList), extent = [0, len(allGenesList[0]), 0, len(allGenesList)], aspect='auto')
    plt.xticks(range(0, len(aminoAcidOrder), 1), rotation=90)
    ax1.set_xticklabels(aminoAcidOrder)
    plt.colorbar(results)

    ###################################################################################     
    #trna given aa
    aminoAcidOrder = get_aa_sort_list(genesDict, typeFlag='aaTrna')
    allGenesList = []
    for x in sortedLociNames:
        try:
            allGenesList.append(sort_amino_acids(allGenesTRNAaa[x], aminoAcidOrder))
        except KeyError:
            pass
    fig = plt.figure()
    ax1 = fig.add_subplot(111)
    results = ax1.matshow(np.array(allGenesList), extent = [0, len(allGenesList[0]), 0, len(allGenesList)], aspect='auto')
    plt.xticks(range(0, len(aminoAcidOrder), 1), rotation=90)
    ax1.set_xticklabels(aminoAcidOrder)
    plt.colorbar(results)
    plt.show()
    
    
    ###################################################################################     
    #writing the data
    ###################################################################################     
    aas = allGenesCodonAA[allGenesCodonAA.keys()[0]].keys()
    outfile = '../Data/codonGivenAA.txt'
    f = open(outfile, 'w')
    for gene in allGenesCodonAA.keys():
        f.write(gene + ',')
        for aa in aas[:-1]:
            f.write(str(allGenesCodonAA[gene][aa]) + ',')
        f.write(str(allGenesCodonAA[gene][aas[-1]]) + '\n')
    f.close()

    aas = allGenesCodonTRNA[allGenesCodonTRNA.keys()[0]].keys()
    outfile = '../Data/codonGivenTRNA.txt'
    f = open(outfile, 'w')
    for gene in allGenesCodonTRNA.keys():
        f.write(gene + ',')
        for aa in aas[:-1]:
            f.write(str(allGenesCodonTRNA[gene][aa]) + ',')
        f.write(str(allGenesCodonTRNA[gene][aas[-1]]) + '\n')
    f.close()
    
    aas = allGenesTRNAaa[allGenesTRNAaa.keys()[0]].keys()
    outfile = '../Data/trnaGivenAA.txt'
    f = open(outfile, 'w')
    for gene in allGenesTRNAaa.keys():
        f.write(gene + ',')
        for aa in aas[:-1]:
            f.write(str(allGenesTRNAaa[gene][aa]) + ',')
        f.write(str(allGenesTRNAaa[gene][aas[-1]]) + '\n')
    f.close()
    
    return 

def compare_correlations():
    genbankFile = '../../CodonAnalysis/data/Genomes/escherichia/escherichia.gb'
    genes, lociNames = pgb.get_coding_sequences(genbankFile, returnLociNames=True)
    #Put them into a dictionary
    genesDict = {}
    for x in range(len(lociNames)):
        genesDict[lociNames[x]] = genes[x]
    
    ###############################################
    #loading some basic dictionaries
    aaDict = saf.get_AA_dict(stop=False)
    codonToAAdict = saf.get_codon_to_AA_dict(stop=False)
    allCodons = codonToAAdict.keys()
    trnaDict = get_trna_dict()
    
    tempTransDict = compileExpressionData.read_protein_expression_data('PNAS-digitalCountsAvg', '../../GeneralFunctions/Data/proteinExpressionP.csv')
    for key in tempTransDict.keys():
        if key not in genesDict.keys():
            del tempTransDict[key]
    tempProteinDict = compileExpressionData.read_protein_expression_data('S-meanProtein', '../../GeneralFunctions/Data/proteinExpressionP.csv')
    for key in tempProteinDict.keys():
        if key not in genesDict.keys():
            del tempProteinDict[key]

    caiWeights = geneMetrics.import_weights_dict('../../GeneralFunctions/Data/geneMetrics/mg1655originalCAIweights.json')

    lengthsDict = {}
    gcDict = {}
    caiDict = {}
    wrightEncDict = {}
    sunEncDict = {}
    infoEncAaDict = {}
    infoEncTrnaDict = {}
    infoEntrnaAaDict = {}
    infoEncAaDictGC = {}
    infoEncTrnaDictGC = {}
    infoEntrnaAaDictGC = {}
    transDict = {}
    proteinDict = {}

    for locusName in genesDict.keys():
        codonList = saf.split_string_into_list_of_triplets(genesDict[locusName])[:-1]
        aaList = translate_codon_list(codonList, codonToAAdict)
        lengthsDict[locusName] = len(codonList)
        gcDict[locusName] = get_gc_percent(genesDict[locusName])
        caiDict[locusName] = calculate_CAI.calculate_CAI(genesDict[locusName], caiWeights)
        temp = enc_wright(codonList)
        if temp <= 61:
            wrightEncDict[locusName] = enc_wright(codonList)
        else:
            wrightEncDict[locusName] = 61
        sunEncDict[locusName] = enc_sun(codonList)
        try:
            transDict[locusName] = tempTransDict[locusName][0]
        except KeyError:
            pass
        try:
            proteinDict[locusName] = tempProteinDict[locusName][0]
        except KeyError:
            pass
        listy = cod_trna_aa_entropy.compute_compressibility_et_al([codonList], gc_model=False, prior_aa=1)
        infoEncAaDict[locusName] = listy[0]['overall_comp']
        infoEncTrnaDict[locusName] = listy[1]['overall_comp']
        infoEntrnaAaDict[locusName] = listy[2]['overall_comp']
        listy = cod_trna_aa_entropy.compute_compressibility_et_al([codonList], gc_model=True, prior_from_gc_model=True)
        infoEncAaDictGC[locusName] = listy[0]['overall_comp']
        infoEncTrnaDictGC[locusName] = listy[1]['overall_comp']
        infoEntrnaAaDictGC[locusName] = listy[2]['overall_comp']
    
    for key in wrightEncDict.keys():
        if wrightEncDict[key] == 'n/a':
            del wrightEncDict[key]

    print stats.spearmanr([caiDict[x] for x in transDict.keys()], [transDict[x] for x in transDict.keys()])
    print stats.spearmanr([caiDict[x] for x in proteinDict.keys()], [proteinDict[x] for x in proteinDict.keys()])

    print 'Transcripts:'
    print stats.spearmanr([sunEncDict[x] for x in transDict.keys()], [transDict[x] for x in transDict.keys()])
    print stats.spearmanr([wrightEncDict[x] for x in transDict.keys() if x in wrightEncDict.keys()], [transDict[x] for x in transDict.keys() if x in wrightEncDict.keys()])

    print stats.spearmanr([infoEncAaDict[x] for x in transDict.keys()], [transDict[x] for x in transDict.keys()])
    print stats.spearmanr([infoEncTrnaDict[x] for x in transDict.keys()], [transDict[x] for x in transDict.keys()])
    print stats.spearmanr([infoEntrnaAaDict[x] for x in transDict.keys()], [transDict[x] for x in transDict.keys()])
    
    print 'Protein:'
    print stats.spearmanr([sunEncDict[x] for x in proteinDict.keys()], [proteinDict[x] for x in proteinDict.keys()])
    print stats.spearmanr([wrightEncDict[x] for x in proteinDict.keys() if x in wrightEncDict.keys()], [proteinDict[x] for x in proteinDict.keys() if x in wrightEncDict.keys()])

    print stats.spearmanr([infoEncAaDict[x] for x in proteinDict.keys()], [proteinDict[x] for x in proteinDict.keys()])
    print stats.spearmanr([infoEncTrnaDict[x] for x in proteinDict.keys()], [proteinDict[x] for x in proteinDict.keys()])
    print stats.spearmanr([infoEntrnaAaDict[x] for x in proteinDict.keys()], [proteinDict[x] for x in proteinDict.keys()])
    
    print 'CAI:'
    print stats.spearmanr([sunEncDict[x] for x in caiDict.keys()], [caiDict[x] for x in caiDict.keys()])
    print stats.spearmanr([wrightEncDict[x] for x in caiDict.keys() if x in wrightEncDict.keys()], [caiDict[x] for x in caiDict.keys() if x in wrightEncDict.keys()])

    print stats.spearmanr([infoEncAaDict[x] for x in caiDict.keys()], [caiDict[x] for x in caiDict.keys()])
    print stats.spearmanr([infoEncTrnaDict[x] for x in caiDict.keys()], [caiDict[x] for x in caiDict.keys()])
    print stats.spearmanr([infoEntrnaAaDict[x] for x in caiDict.keys()], [caiDict[x] for x in caiDict.keys()])
    
    plt.close()
    
    plt.figure()
    plt.scatter([wrightEncDict[x] for x in wrightEncDict.keys()], [infoEncAaDict[x] for x in wrightEncDict.keys()])
    plt.xlabel('Wright\'s ENC')
    plt.ylabel('Info Bias')
    
    plt.figure()
    plt.scatter([lengthsDict[x] for x in wrightEncDict.keys()], [wrightEncDict[x] for x in wrightEncDict.keys()])
    plt.xlabel('Length in codons')
    plt.ylabel('Wright\'s ENC')
    
    plt.figure()
    plt.scatter([lengthsDict[x] for x in infoEncAaDict.keys()], [infoEncAaDict[x] for x in infoEncAaDict.keys()])
    plt.xlabel('Length in codons')
    plt.ylabel('Info Bias')
    
    plt.figure()
    plt.scatter([infoEncTrnaDict[x] for x in infoEncAaDict.keys()], [infoEntrnaAaDict[x] for x in infoEncAaDict.keys()], c=[infoEncAaDict[x] for x in infoEncAaDict.keys()])
    plt.xlabel('Info bias: Codon given tRNA')
    plt.ylabel('Info Bias: tRNA given AA')
    cb = plt.colorbar()
    cb.set_label('Info Bias: Codon given AA')
    
    plt.figure()
    plt.scatter([infoEncTrnaDict[x] for x in transDict.keys()], [infoEntrnaAaDict[x] for x in transDict.keys()], c=[log(transDict[x]) for x in transDict.keys()])
    plt.xlabel('Info bias: Codon given tRNA')
    plt.ylabel('Info Bias: tRNA given AA')
    cb = plt.colorbar()
    cb.set_label('Log transcript abundance')
    
    plt.figure()
    plt.scatter([infoEncAaDict[x] for x in infoEncAaDict.keys()], [infoEntrnaAaDictGC[x] for x in infoEncAaDict.keys()], c=[gcDict[x] for x in infoEncAaDict.keys()])
    plt.xlabel('Info bias: Non-GC model')
    plt.ylabel('Info Bias: GC model')
    cb = plt.colorbar()
    cb.set_label('GC percent')
    
    plt.figure()
    plt.scatter([gcDict[x] for x in wrightEncDict.keys()], [wrightEncDict[x] for x in wrightEncDict.keys()])
    plt.xlabel('GC percent')
    plt.ylabel('Wright\'s Enc')
    
    plt.figure()
    plt.scatter([gcDict[x] for x in wrightEncDict.keys()], [infoEncAaDict[x] for x in wrightEncDict.keys()])
    plt.xlabel('GC percent')
    plt.ylabel('Info Bias: Non-GC model')
    
    plt.figure()
    plt.scatter([gcDict[x] for x in wrightEncDict.keys()], [infoEncAaDictGC[x] for x in wrightEncDict.keys()])
    plt.xlabel('GC percent')
    plt.ylabel('Info Bias: GC model')
    
    plt.show()

    return 

def pca_loadings_vs_expression(CAIfile, PCAloadingsFile, CompressFile):
    dataDictTranscript = compileExpressionData.read_protein_expression_data('PNAS-digitalCountsAvg', '../../GeneralFunctions/Data/proteinExpressionP.csv')
    dataDictProtein = compileExpressionData.read_protein_expression_data('S-meanProtein', '../../GeneralFunctions/Data/proteinExpressionP.csv')
    dataDictCAI = {}
    
    f=open(CAIfile, 'r').readlines()
    for x in f:
        line = x.split(',')
        dataDictCAI[line[0]] = float(line[1])
    
    tempList = []
    f=open(PCAloadingsFile, 'r').readlines()
    for x in f:
        line = x.split(' ')
        tempList.append(float(line[0]))
    
    dataDictPCA = {}
    f=open(CompressFile, 'r').readlines()
    for x in range(len(f)):
        line = f[x].split(',')
        dataDictPCA[line[0]] = tempList[x]
    
    dataTrans = [dataDictTranscript[x][0] for x in dataDictTranscript.keys() if x in dataDictPCA.keys()]
    dataTransPCA = [dataDictPCA[x] for x in dataDictTranscript.keys() if x in dataDictPCA.keys()]
    print stats.spearmanr(dataTrans, dataTransPCA)
    dataPro = [dataDictProtein[x][0] for x in dataDictProtein.keys() if x in dataDictPCA.keys()]
    dataProPCA = [dataDictPCA[x] for x in dataDictProtein.keys() if x in dataDictPCA.keys()]
    print stats.spearmanr(dataPro, dataProPCA)
    dataCAI = [dataDictCAI[x] for x in dataDictCAI.keys() if x in dataDictPCA.keys()]
    dataPCA = [dataDictPCA[x] for x in dataDictCAI.keys() if x in dataDictPCA.keys()]
    print stats.spearmanr(dataCAI, dataPCA)
    return

def column_wise_methods():
    genbankFile = '../../CodonAnalysis/data/Genomes/escherichia/escherichia.gb'
#    genes, lociNames = pgb.get_coding_sequences(genbankFile, returnLociNames=True)
    
    refSetFile = '../../CodonAnalysis/data/referenceSets/top25PercentProtein'
#    refSetFile = '../../CodonAnalysis/data/referenceSets/escherichiaOriginalCAIset'
    genes = pgb.get_coding_sequences(genbankFile, subsetFileLocation=refSetFile)


    codonListOfLists = [saf.split_string_into_list_of_triplets(gene)[:-1] for gene in genes]
    columnListOfLists = [[] for x in range(300)]
    for x in range(300):
        for gene in codonListOfLists:
            try:
                columnListOfLists[x].append(gene[x])
            except IndexError:
                pass
    aaCodonList = []
    aaTrnaList = []
    trnaCodonList = []
    originalEncList = []
    for codonList in columnListOfLists:
        listy = cod_trna_aa_entropy.compute_compressibility_et_al([codonList], gc_model=True, prior_from_gc_model=True)
        aaCodonList.append(listy[0]['overall_comp'])
        trnaCodonList.append(listy[1]['overall_comp'])
        aaTrnaList.append(listy[2]['overall_comp'])
        originalEncList.append(enc_wright(codonList))
    fig1=plt.figure(1)
    ax1=fig1.add_subplot(111)
    ax1.plot(aaCodonList[1:])
    
    fig2=plt.figure(2)
    ax2=fig2.add_subplot(111)
    ax2.plot(trnaCodonList[1:])
    
    fig3=plt.figure(3)
    ax3=fig3.add_subplot(111)
    ax3.plot(aaTrnaList[1:])
    
    fig4=plt.figure(4)
    ax4=fig4.add_subplot(111)
    ax4.plot(originalEncList[1:])
    
    refSetFile = '../../CodonAnalysis/data/referenceSets/bottom25PercentProtein'
    genes = pgb.get_coding_sequences(genbankFile, subsetFileLocation=refSetFile)


    codonListOfLists = [saf.split_string_into_list_of_triplets(gene)[:-1] for gene in genes]
    columnListOfLists = [[] for x in range(300)]
    for x in range(300):
        for gene in codonListOfLists:
            try:
                columnListOfLists[x].append(gene[x])
            except IndexError:
                pass
    aaCodonList = []
    aaTrnaList = []
    trnaCodonList = []
    originalEncList = []
    for codonList in columnListOfLists:
        listy = cod_trna_aa_entropy.compute_compressibility_et_al([codonList], gc_model=True, prior_from_gc_model=True)
        aaCodonList.append(listy[0]['overall_comp'])
        trnaCodonList.append(listy[1]['overall_comp'])
        aaTrnaList.append(listy[2]['overall_comp'])
        originalEncList.append(enc_wright(codonList))
    
    ax1.plot(aaCodonList[1:])
    ax2.plot(trnaCodonList[1:])
    ax3.plot(aaTrnaList[1:])
    ax4.plot(originalEncList[1:])
    
    
    
    plt.show()
    return


########################################################################################
#very basic support functions
########################################################################################
def split_tuple_dict(aaDict, geneAAdict, metric='compress'):
    if metric == 'compress':
        tupleValue = 0
    else:
        tupleValue = 1
    copyDict = {}
    for x in geneAAdict.keys():
        copyDict[x] = geneAAdict[x][tupleValue]
    return copyDict

def sort_amino_acids(dictionary, aaList):
    listy = []
    for aa in aaList:
        listy.append(dictionary[aa])
    return listy

def get_aa_sort_list(genesDict, typeFlag):
    codonList = [saf.split_string_into_list_of_triplets(genesDict[gene])[:-1] for gene in genesDict.keys()]
    codonList = [inner for outer in codonList for inner in outer]
#    listy = cod_trna_aa_entropy.compute_compressibility_et_al([codonList], gc_model=True, prior_from_gc_model=True)
    listy = cod_trna_aa_entropy.compute_compressibility_et_al([codonList], gc_model=False, prior_aa=1)
    aaCodonDict = listy[0]
    trnaCodonDict = listy[1]
    aaTRNAdict = listy[2]
    if typeFlag == 'aaCodon':
        sortedData = sorted(aaCodonDict['x_to_comp_effno'].iteritems(), key=operator.itemgetter(1))
    elif typeFlag == 'aaTrna':
        sortedData = sorted(aaTRNAdict['x_to_comp_effno'].iteritems(), key=operator.itemgetter(1))
    elif typeFlag == 'trnaCodon':
        sortedData = sorted(trnaCodonDict['x_to_comp_effno'].iteritems(), key=operator.itemgetter(1))

    aaList = [x[0] for x in sortedData]
    return aaList

def translate_codon_list(codonList, codonToAAdict):
    aminoAcidList = [codonToAAdict[x] for x in codonList]
    return aminoAcidList

def get_trna_dict():
    trnaDict = {
            'Phe1': ['TTT', 'TTC'],
            'Leu1': ['TTA', 'TTG'],
            'Leu2': ['CTT', 'CTC'],
            'Leu3': ['CTA', 'CTG'],
            'Ile1': ['ATT', 'ATC'],
            'Ile2': ['ATA'],
            'Met1': ['ATG'],
            'Val1': ['GTT', 'GTC'],
            'Val2': ['GTA', 'GTG'],
            'Ser1': ['TCT', 'TCC'],
            'Ser2': ['TCA', 'TCG'],
            'Ser3': ['AGT', 'AGC'],
            'Pro1': ['CCT', 'CCC'],
            'Pro2': ['CCA', 'CCG'],
            'Thr1': ['ACT', 'ACC'],
            'Thr2': ['ACA', 'ACG'],
            'Ala1': ['GCT', 'GCC'], 
            'Ala2': ['GCA', 'GCG'],
            'Tyr1': ['TAT', 'TAC'],
            'His1': ['CAT', 'CAC'],
            'Glu1': ['CAA', 'CAG'],
            'Asp1': ['AAT', 'AAC'],
            'Lys1': ['AAA', 'AAG'],
            'AsA1': ['GAT', 'GAC'],
            'GlA1': ['GAA', 'GAG'],
            'Cys1': ['TGT', 'TGC'],
            'Trp1': ['TGG'],
            'Arg1': ['CGT', 'CGC'],
            'Arg2': ['CGA', 'CGG'],
            'Arg3': ['AGA', 'AGG'],
            'Gly1': ['GGT', 'GGC'],
            'Gly2': ['GGA', 'GGG']
            }
    return trnaDict

def get_codon_to_trna_dict(trnaDict):
    dicty = {}
    for trna in trnaDict.keys():
        for codon in trnaDict[trna]:
            dicty[codon] = trna
    return dicty

def get_aa_to_trna_dict(aaToCodonDict, codonToTRNAdict):
    dicty = {}
    for aa in aaToCodonDict.keys():
        listy = []
        for codon in aaToCodonDict[aa]:
            listy.append(codonToTRNAdict[codon])
        dicty[aa] = list(set(listy))
    return dicty




########################################################################################
########################################################################################
#deprecated/unused functions
########################################################################################
########################################################################################

########################################################################################
#Generating and testing basic strings
def generate_test_sequence(aaProbabilityDict, numberToRemove, length):
    aaDict = saf.get_AA_dict(stop=False)
    for x in range(numberToRemove):
        changed = False
        while not changed:
            aa = random.choice(aaDict.keys())
            if len(aaDict[aa]) > 1:
                aaDict[aa].pop()
                changed=True
    codonList = []
    aminoAcidList = []
    for aminoAcid in aaDict.keys():
        if int(length*aaProbabilityDict[aminoAcid]) > 0:
            for codonChoice in range(int(length*aaProbabilityDict[aminoAcid])):
                codonList.append(random.choice(aaDict[aminoAcid]))
                aminoAcidList.append(aminoAcid)
    return codonList, aminoAcidList

def basic_tests():
    genbankFile = '../../CodonAnalysis/data/Genomes/escherichia/escherichia.gb'
    genes, lociNames = pgb.get_coding_sequences(genbankFile, returnLociNames=True)
    codonToAAdict = saf.get_codon_to_AA_dict(stop=False)
    aaListOfLists = [translate_codon_list(saf.split_string_into_list_of_triplets(x)[:-1], codonToAAdict) for x in genes]

    numberToRemove = 10 
    iterations = 100
    aaDict = saf.get_AA_dict(stop=False)
    allCodons = [inner for outer in aaDict.values() for inner in outer] 
    #uniform amino acid structure###################################
    aaProbabilityDict = saf.get_AA_dict(stop=False)
    for aa in aaProbabilityDict.keys():
        aaProbabilityDict[aa] = 1/float(20)
    actualGeneAAfreq = random.choice(aaListOfLists)
    for aa in aaProbabilityDict.keys():
        aaProbabilityDict[aa] = actualGeneAAfreq.count(aa)/float(len(actualGeneAAfreq))
    
    print aaProbabilityDict

    lengthRange = range(100, 2100, 100)
    prior = 0.0001
    
    encWright = {}
    encSun = {}
    encSunInfo = {}
    ourInfoCompress = {}
    ourInfoUnweighted = {}
    for lengthToTest in lengthRange:
        encWright[lengthToTest] = []
        encSun[lengthToTest] = []
        encSunInfo[lengthToTest] = []
        ourInfoCompress[lengthToTest] = []
        ourInfoUnweighted[lengthToTest] = []
        
        for x in range(iterations):
            codonList, aaList = generate_test_sequence(aaProbabilityDict, numberToRemove, lengthToTest)
            infoTheoryCompress, infoTheoryUnweighted = codon_entropy.compute_entropies([codonList], [aaList], allCodons, aaDict, 0, prior)
            encWright[lengthToTest].append(enc_wright(codonList))
            encSun[lengthToTest].append(enc_sun(codonList, prior=prior))
            encSunInfo[lengthToTest].append(enc_sun(codonList, publishedFcalculation=False, publishedWeightScheme=False, prior=prior))
            ourInfoCompress[lengthToTest].append(infoTheoryCompress)
            ourInfoUnweighted[lengthToTest].append(infoTheoryUnweighted)
    
    plt.close()
    plt.figure()
    xaxis=[]
    yaxis= []
    err = []
    for lengthToTest in lengthRange:
        xaxis.append(lengthToTest)
        yaxis.append(np.mean(encWright[lengthToTest]))
        err.append(np.std(encWright[lengthToTest]))
    plt.errorbar(xaxis, yaxis, err)
    plt.axhline(61-numberToRemove, linewidth=3, color='k')
    plt.xlim(0, lengthRange[-1] + 100)
    ylimit = plt.ylim()
   
    plt.figure()
    xaxis=[]
    yaxis= []
    err = []
    for lengthToTest in lengthRange:
        xaxis.append(lengthToTest)
        yaxis.append(np.mean(encSun[lengthToTest]))
        err.append(np.std(encSun[lengthToTest]))
    plt.errorbar(xaxis, yaxis, err, color='g')
    xaxis=[]
    yaxis= []
    err = []
    for lengthToTest in lengthRange:
        xaxis.append(lengthToTest)
        yaxis.append(np.mean(encSunInfo[lengthToTest]))
        err.append(np.std(encSunInfo[lengthToTest]))
    plt.errorbar(xaxis, yaxis, err, color='r')
    plt.axhline(61-numberToRemove, linewidth=3, color='k')
    plt.ylim(ylimit)

    plt.figure()
    xaxis=[]
    yaxis= []
    err = []
    for lengthToTest in lengthRange:
        xaxis.append(lengthToTest)
        yaxis.append(np.mean(ourInfoCompress[lengthToTest]))
        err.append(np.std(ourInfoCompress[lengthToTest]))
    plt.errorbar(xaxis, yaxis, err, color='c')
    xaxis=[]
    yaxis= []
    err = []
    for lengthToTest in lengthRange:
        xaxis.append(lengthToTest)
        yaxis.append(np.mean(ourInfoUnweighted[lengthToTest]))
        err.append(np.std(ourInfoUnweighted[lengthToTest]))
    plt.errorbar(xaxis, yaxis, err, color='m')
    plt.axhline(61-numberToRemove, linewidth=3, color='k')
    plt.ylim(ylimit)
    
    plt.figure()
    plt.hist(encWright[1000], label='Wright ' + str(round(np.mean(encWright[1000]),2)), alpha=0.7)
    plt.hist(encSun[1000], label='Sun ' + str(round(np.mean(encSun[1000]), 2)), alpha=0.7)
    plt.hist(encSunInfo[1000], label='SunInfo ' + str(np.mean(encSunInfo[1000])))
    plt.hist(ourInfoCompress[1000], label='InfoCompress ' + str(round(np.mean(ourInfoCompress[1000]),2)), alpha=0.7)
    plt.hist(ourInfoUnweighted[1000], label='InfoUnweight ' + str(round(np.mean(ourInfoUnweighted[1000]),2)), alpha=0.7)
    plt.axvline(61-numberToRemove, linewidth=3, color='k')
    plt.legend(loc=2)
    plt.show()
    return

def simple_test():
    aaDict = saf.get_AA_dict(stop=False)
#    listy = aaDict['serine'] * 10 + aaDict['serine'][:5] * 10
    listy = aaDict['serine'][:1] * 30 + aaDict['serine'][:5] * 10
    possibilities = {}
    for x in listy:
        try:
            possibilities[x] += 1
        except:
            possibilities[x] = 1
    probabilities = [possibilities[x]/float(sum(possibilities.values())) for x in possibilities.keys()]
    print possibilities
    infoTheory = [-1.* x * log(x,2) for x in probabilities]
    print 2**sum(infoTheory)
    print 1/sum([x**2 for x in probabilities])
    return



########################################################################################
#getting distance csvs to make a similarity network

def get_euclidean_distance(dicty1, dicty2):
    assert dicty1.keys() == dicty2.keys(), "The keys of the dictionaries are different!"
    dist = np.linalg.norm(np.array(dicty1.values())-np.array(dicty2.values()))
    return dist

def get_cosine_distance(dicty1, dicty2):
    assert dicty1.keys() == dicty2.keys(), "The keys of the dictionaries are different!"
    dist = distance.cosine(np.array(dicty1.values()), np.array(dicty2.values()))
    return dist

def get_distance_csvs():
    genbankFile = '../../CodonAnalysis/data/Genomes/escherichia/escherichia.gb'
    genes, lociNames = pgb.get_coding_sequences(genbankFile, returnLociNames=True)
    dataDict = {}
    genome = {}
    for x in range(len(lociNames)):
        genome[lociNames[x]] = genes[x]
    aaDict = saf.get_AA_dict(stop=False)
    codonToAAdict = saf.get_codon_to_AA_dict(stop=False)
    allCodons = aaDict.values()
    allCodons = [inner for outer in allCodons for inner in outer] 
    for gene in genome.keys():
        codonList = saf.split_string_into_list_of_triplets(genome[gene])[1:-1]
        aaList = translate_codon_list(codonList, codonToAAdict)
        geneAAdict = codon_entropy.compute_entropies([codonList], [aaList], allCodons, aaDict, 0, 1, aa_dict_flag=True)
        dataDict[gene] = geneAAdict
#    fEuc = open('../Data/euclideanDistance.csv', 'w')
#    fCos = open('../Data/cosineSimilarity.csv', 'w')
#    fDict = open('../Data/aaEncVectors.csv', 'w')
    for aa in dataDict[dataDict.keys()[0]]:
        print aa
#    for gene1 in dataDict.keys():
#        fDict.write(gene1 + ',')
#        for aa in dataDict[gene1]:
#            if aa == 'tyrosine':
#                fDict.write(str(dataDict[gene1][aa]) + '\n')
#            else:
#                fDict.write(str(dataDict[gene1][aa]) + ',')
#        for gene2 in dataDict.keys():
#            if gene2 != gene1:
#                euclideanDistance = get_euclidean_distance(dataDict[gene1], dataDict[gene2])
#                cosineSimilarity = get_cosine_distance(dataDict[gene1], dataDict[gene2])
#                fEuc.write(gene1 + ',' + gene2 + ',' + str(euclideanDistance) + '\n')
#                fCos.write(gene1 + ',' + gene2 + ',' + str(cosineSimilarity) + '\n')
#        del dataDict[gene1]
    return
